Exchange-Correlation Energy Functional Based on the Airy-Gas Reference System 
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In recent work, generalized gradient approximations (GGA's) have been constructed from the 
energy density of the Airy gas for exchange but not for correlation. We report the random phase 
approximation (RPA) conventional correlation energy density of the Airy gas, the simplest edge 
electron gas, in which the auxiliary noninteracting electrons experience a linear potential. By 
fitting the Airy-gas RPA exchange-correlation energy density and making an accurate short-range 
correction to RPA, we propose a simple beyond-RPA GGA density functional ("ARPA-I-") for 
the exchange-correlation energy. Our functional, tested for jellium surfaces, atoms, molecules and 
solids, improves mildly over the local spin density approximation for atomization energies and lattice 
constants without much worsening the already-good surface exchange-correlation energies. 

PACS numbers: 71.10.Ca,71.15.Mb,71.45.Gm 



I. INTRODUCTION 

In Kohn-Sham density functional theory^, the ground- 
state density and energy of interacting electrons in a 
scalar external potential v{y) are computed efficiently 
via a selfconsistent calculation for an auxiliary system 
of noninteracting electrons in a scalar effective potential 
Vef / (r) . Once the exchange-correlation energy as a func- 
tional of the electron density has been approximated, its 
functional derivative provides the exchange-correlation 
contribution to v^f / (r) . By itself, the deviation of Ue/ / (r) 
from the constant chemical potential determines the elec- 
tron density and thus the correlation energy. Typical 
approximations are designed to be exact for a reference 
system, most often the uniform electron gas in which 
the auxiliary noninteracting electrons see a constant or 
uniform v^f f . Sometimes additional exact constraints or 
fits to experiment are also built into the approximation. 
Recently Kohn and Mattsson^ have proposed as a more 
realistic reference system the edge electron gas, in which 
Vef f (r) varies more or less linearly near the edge surface 
of the density. While the uniform gas could be (and is) 
a good reference for a bulk solid, the edge electron gas 
could be at least as good for a bulk solid and better for 
solid surfaces, molecules, and atoms, which have regions 
where the electron density evanesces. 

The edge surface of any electron system is defined^ 
by Veff{r) = /i, where Ue//(r) is the exact Kohn-Sham^ 
(KS) effective potential and fi is the chemical potential. 
Outside this classical turning surface, all noninteracting 
electrons tunnel into a barrier. The simplest example of 
an edge electron gas is the Airy gas, where any electron 
feels a linear effective potential and thus the normalized 
one-particle eigenfunctions are proportional to the Airy 
function. The Airy gas has not only a surface-like region, 
but also a region of high and slowly-varying (Thomas- 
Fermi-like) electron density where the local density ap- 
proximation (with uniform-gas input) is accurate^i^ for 
the noninteracting kinetic, exchange, and correlation en- 
ergy densities. 

The Airy gas has appeared before in density functional 



theory: (1) The effective finite-linear-potential model 
gives remarkably good results for the jellium surface 
problem, where the orbitals of this model are approxi- 
mated with plane waves inside the bulk, Airy functions 
near the surface, and exponential functions far in the 
vacuum^ (2) Baltinl constructed a generalized gradi- 
ent approximation (GGA) for the orbital kinetic energy 
from the Airy-gas kinetic energy density, but his approx- 
imation does not recover the second-order gradient ex- 
pansion for the kinetic energy density^i^ and is poor for 
atoms and molecule o^°i^^ . However, the kinetic energy 
density of the Airy gag-'^^ can still be a starting point for 
construction of GGA kinetic energy functionals that can 
be more accurate for atoms, molecules, jellium clusters, 
and jellium surfaceaii^i^. The trick is to fit a GGA plus 
a V^n term integrating to zero to the Airy-gas kinetic 
energy density. 

The exchange energy density of the Airy gas^ was 
fittedi^ii^ with a function dependent on the density and 
its gradient. Thus, Vitos et. alM developed a GGA ex- 
change energy functional (LAG or local Airy-gas GGA) 
that was used with the local spin-density approximation 
(LSDA) correlation energy. This exchange-correlation 
(xc) energy functional gives results for atoms very close 
to, but better than, the LSDA ones, and its accuracy 
for atomization energy of diatomic molecules is similar 
to that of the PBE GGAi^, while for bulk systems the 
resuhs of LAG GGA are close to the PBEsol GGA^^ and 
to experimental values. However, the jellium xc surface 
energies of LAG are far too low (lower even than those of 
the PBE GGA). Armiento and Mattssoi*i^iii proposed 
an xc energy functional (AM05 GGA) using a better fit 
for the Airy gas exchange energy density and a correla- 
tion energy functional constructed such that the AM05 
xc jellium surface energies fit the RPA-l-iS, values (RPA 
plus a GGA short-range correction). AMOS is also based 
on the subsystem functional approach^-., which permits 
an interpolation between a uniform-gas reference for the 
bulk of a solid and an Airy-gas reference for the sur- 
face. (Since the Airy-gas reference system by itself pro- 
vides such an interpolation, we make no further inter- 
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polation here.) AMOS slightly improves the accuracy of 
LAG GGA for bulk systems. 

Because the correlation energy density of the Airy gas 
was unknown, the LAG GGA and AMOS GGA used 
in their construction only the Airy-gas exchange energy 
density. In this paper we compute the correlation energy 
density of the Airy gas in the random phase approxima- 
tion (RPA), and fit it to a GGA (ARPA). As in Refs.^ 
and^**, our fit is made without regard to exact constraints 
on Exc[n^,ni]. The Airy gas is a system of delocalized 
electrons where the self-interaction correction has no ef- 
fect, and where the GGA correction^ to the integrated 
RPA energy should be accurate. Our functional, includ- 
ing this GGA correction to RPA, will be called ARPA+. 

Unlike energies, energy densities of non-uniform sys- 
tems are not unique. It is not clear to us that the conven- 
tional choice for the exchange-correlation energy density 
(made in RefsJ^ii^, and here) is optimal. It is not our 
intention here to either endorse or criticize this choice, 
but simply to see what GGA is obtained from the Airy- 
gas reference system within a consistent implementation 
for correlation as well as exchange. 

AMOS, PBEsol, and ARPA -I- are of special interest as 
candidates for a " GGA for solids" providing better lattice 
constants and surface energies than standard GGA's like 
PBE, possibly at the cost of a worsened description of 
atoms and molecules. There have been several recent 
articles commenting on or testing for solids the LAG, 
AMOS, and PBEsol GGA's^i2i^i23^^ 

Our paper is organized as follows. In section|TTl we pro- 
pose a simple model for the Airy gas. In section ITTTl we 
construct the ARPA -I- GGA xc energy functional from 
our Airy gas model. In section IIVI we test the ARPA+ 
GGA for atoms, molecules, jellium surfaces and bulk 
solids. In section |Vl we summarize our conclusions. 



the equation 

with the boundary conditions 

0j(-c^) =0,(L) = 0. 
They are given by the Airy functions 



0,(z) = aAi(---^) 



(3) 



(4) 



(5) 



where e = (F^/2)^/'^ is the Airy gas characteristic en- 
ergy scale, a is the normalization constant, and ej is the 
j-th eigenvalue calculated from the boundary condition 
4>j{L) = 0. The Airy gas density is 



l(z)=^0,2(^)|^^.|/^_ 



(6) 



We recall that all 3D states with energy up to = are 
occupied. Thus the Airy gas is completely determined 
by the length I and the energy e. 

In the limit L/l oo, the normalization constant is^ 



rl/2 



and the eigenvalues are^ 



So, the density of the Airy gas is 

n{z) = r^no{f]), T] = z/l, 



(7) 



(8) 



(9) 



II. THE AIRY GAS MODEL 

The simplest example of an edge electron gas is the 
Airy gas that is translationally invariant in the plane of 
the surface {z = 0) and has the effective potential3.i2S. 



where 



-Fz, -oo < z < L (F > 0) 
oo, z > L {L/l — > oo). 



(1) 



Here F ~ \dVf,j f{z) / dz\ is the slope of the effective po- 
tential and the characteristic length scale 



/ = (2F)-i/3 



(2) 



is approximately the edge region thickness^. (Unless oth- 
erwise stated, atomic units are used throughout, i.e., 
= h = rrie = 1.) 

The KS orbitals are ^'j,k||(r) = ^i'j (z)-^e**'iirii , where 
k|| and r|| are the wavevector and the position vector par- 
allel to the plane of the surface, A is the cross-sectional 
area, and the orthonormal eigenfunctions (l)j{z) satisfy 
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(10) 



Let us consider a model for the Airy gas that is 
described by Eqs. (U) - ([6]), but instead of choosing 
L// — » oo we take L/l = 20 for computational conve- 
nience. Such a system has 19 occupied orbitals 0j(z) 
and can accurately describe the Airy gas. The normal- 
ization constants of Eq. ([S]) and the eigenvalues ej are 
computed numerically. Such an approach is similar to 
jellium slabs that are described by a finite number of 
occupied orbitals in the z-direction and that can accu- 
rately predict the surface energies of semi-infinite jellium 
surfaces- 

We select three values F — 0.1, F — O.S, and F = 1 
for the slope of the effective potential. The accuracy of 
the model does not depend on the F value. In Fig. [T]we 
show the densities of the Airy gas and of our Airy gas 
model for the chosen values of the slope F. We see the 
exact Airy gas densities and the modeled ones can not 
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FIG. 1: Electron density (electrons/bohr^) of the Airy gas 
and of our model versus z (bohr), for several slopes of the 
effective potential (F — 0.1 making / — 1.710, F — 0.5 making 
/ = 1.000, F = 1 making I = 0.793). The edge is at 2 = 0. 



FIG. 2: Reduced gradient s{z) versus z, of the Airy gas 
and our model, for several slopes of the effective potential 
(F = 0.1, 0.5, and 1.) The edge is at 2 = 0. 



be distinguished until z ^ L = 20 ■ I where the densities 
of our model have oscillations until they vanish. 

Important ingredients of any GGA functional are the 
density n(r) and the reduced density gradient 



3(r) = |Vn(r)|/[2fcF(r)n(r)] 



(11) 



where fcF(r) = (37r^ri(r))^/'^ is the Fermi wavevector. 
(The dimensionless density gradient s(r) measures the 
variation of the density over a Fermi wavelength Xp — 
2TT/kF.) In Fig. [2] we compare the reduced gradients of 
our model and of the exact Airy gas. Up to s = 2, the 
model nicely matches the exact Airy gas, and it is ac- 
curate for any value of s. (We note that s values bigger 
than 3 are found in the tail of an atom or molecule, where 
the electron density is negligible. We also note that in 
most bulk solids the maximum^^ value of the reduced 
gradient is smaller than 2.) Figs. [T] and [2] demonstrate 
that our model is accurate, and thus we can use it for the 
calculation of the Airy gas correlation energy. 



III. RPA CORRELATION ENERGY DENSITY 
OF THE AIRY GAS, AND THE CONSTRUCTION 
OF THE ARPA+ GGA 

The conventional xc energy density at a point is ne^cy 
where n is the local electron density and e^c is the con- 
ventional xc energy per particle. Let us consider the spin- 
unpolarized Airy gas model with the edge plane at z = 0. 
Using its translational invariance in a plane perpendicu- 
lar to the z axis, and the so-called adiabatic-connection 
fluctuation-dissipation theoieu^'^^^ (ACFDT), the 
exact expression for the conventional xc energy per par- 



ticle at point z is 



26.27.28 



(2^)2 



dz u(z,i,g||)[- 



Tm{z) 



dX 



dujx^{z,z;q\\,iuj) - S(z - z)], (12) 



where q|| is the wavevector parallel to the surface, and 
and V are the two-dimensional Fourier transforms of 
the interacting density response function at the coupling 
strength A and of the Coulomb potential respectively. 
The substitution of with the non-interacting density 
response function into Eg. p^ yields the exact ex{z) 
(expressible in terms of occupied orbitals only, although 
x'^ requires also the unoccupied orbitals). The density 
response function obeys the screening integral Dyson-like 
equation"^ 



X°(r, r',uj) + y dridr2X°(r, ri,w) 



xK(ri,r2) + /4[n](ri,r2,w)}x^(r2,r',c<.), (13) 



where w (ri,r2) = A/|ri 
Sv^^[n](ri,u)/Sn{r2,uj) is 



- ral and /^Jn] (ri, ra, w) = 
the exact xc kernel. Here 
v^^[n\ is the exact frequency-dependent xc potential at 
coupling strength A. Obviously, the exact xc kernel is un- 
known and it has to be approximated. Approximations 
of the xc kernel are usually constructed from the uniform 
electron gas^^^^^, and have not been tested sufRciently 
for nonuniform systems. When /^^[n](r, r'; w) is taken 
to be zero, Eq. reduces to the RPA. The RPA xc 
hole density is exact at large interelectronic separations 
such that it can correctly describe the xc hole density of 
an electron far outside of a jellium surface^**, and its on- 
top hole is finite and well described by the LSDA-RPA3 
on-top hole in the case of a jellium surface^. 

Eqs. and can be generalized^^ for systems 

with any relative spin polarization 



c 



(14) 



where and rij are the spin densities, n| +ni = n. Thus 
for the Airy gas model, we choose to calculate the RPA 
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correlation energy per particle at point z, from Eqs. p2 
and (|13p . and to add the RPA+ short-range correction: 



RPA+ _ j^RPA 



GGA 

xc 



GGA~RPA 



(15) 



where E^^"^ is the PBE GGAi^ xc energy, and 
E'^^'^-^P^ is the PBE-RPA GGA xc energyi^. The ex- 
change contribution and the long-range correlation con- 
tribution cancel out of the bracketed term in Eq. (15), 
leaving only short-range correlation. Because the self- 
interaction correction is not important for the Airy gas, 
Eq. (|15p will give nearly the exact correlation energy of 
the Airy gas. 

For the numerical evaluation of Eqs. and 
we follow the method described in Refsi^ and2^, but in- 
stead of using the double- and single-cosine representa- 
tions of the density response function and the density re- 
spectively, we use a grid on the z-axis for x^(-2^) -2j 9|h *'^) 
and n{z). We find that the first 50 unoccupied orbitals 
<j)j{z) are enough for an accurate calculation. (Our grid 
on the z-axis can accurately describe the occupied and 
the first 50 unoccupied orbitals^). 

The exchange energy for a spin-polarized system may 
be evaluated from the spin-unpolarized version using the 
spin-scaling relation^^: 

E^[n^,n^] = i{^,[2nt]+^,[2nj}, (16) 

and thus we only need to consider the spin-unpolarized 
case. We fit the exchange energy per particle of the Airy 
gas model, using the non-linear least-square Levenberg- 
Marquardt method'^^, with the following expression 

e^{n{v))^e^J^^{n{v))F^{s{v)), (17) 

where e^^^^ = —Skp/iir and the enhancement factor is 



ais" 



(1 + ass"^)" 



1 - 055°'' + ays" 
1 -I- ags'^io 



(18) 



where ai = 0.041106, 03 = 2.626712, as = 0.092070, 
a4 = 0.657946 are the parameters found in Ref.^'^, and 
a5 = 133.983631, ag = 3.217063, 07 = 136.707378, ag = 
3.223476, ag = 2.675484, aw = 3.473804 are parameters 
found from our fitting procedure. Eq. (jl7p recovers the 
correct LSDA for the uniform electron gas, and fits well 
the Airy gas exchange energy per particle for s < 20. 
( Values of s bigger than 20 are found only when the 
density is negligible. We recall that LAA of Reff^^ is a 
better fit than LAG or far outside the edge.) 

In Fig. [3] we show {e^ — e^^^^)/ex versus the re- 
duced gradient s for several approximations. The Airy 
gas curve, as well as our Airy gas model curve, have a 
negative region around s ~ 0.5 that was not taken into 
account by the LAG GGA and AM05 GGA. We find 



this fine feature only because we plot {e^ 



^LSDA 



instead of ex- (This feature can also be seen in the inset 
of Fig. 1 of Ref.^^, but it was not taken into account 
in the construction of AM05.) The second term of the 
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FIG. 3: (e^ — )/ex versus the reduced gradient s for the 

Airy gas model, the Airy gas, the LAG GGAf^, and for our 
fit (see Eqs. ([17]) and d])). The "Airy gas" curve uses 
and e^^^^ of the Airy gas, whereas the other curves use 
and ei/^^-* of our model for the Airy gas. The AMOS GGA^4, 
not shown in the figure, has the same behavior as the LAG 
GGA. 



right-hand-side of Eq. ^TE\i models the exact behavior 
at small reduced gradients, whereas the first term of the 
right-hand-side of Eq. has the same form as the 

parametrization proposed in Refii^. We observe that our 
fit (Eqs. (|17p and ^TE\i ) is very close to the exact Airy gas 
model as well as to the exact Airy gas exchange energy 
per particle. 

We fit the RPA correlation energy per particle of 
the Airy gas of any spin polarization with the follow- 
ing expression, using again the non-linear least-square 
Levenberg-Marquardt method'^^ 

6^«^^(r„ C, .e) = e^'^^-^^^{r,,OFdsc), (19) 

where is the local Wigner-Seitz radius [n — 3/(47rrf ) — 
fc|./37r^], C is the relative spin polarization of Eq. (|14p . 
^LSDA-RPA ^Yie RPA correlation energy per particle 

and 



of the uniform electron gas (see Ref4£) 

Sc(r)-^|Vn(r)|/[2(3^2)i/3„(j.)7.9/6]^ 



(20) 



with (j) = [{1 + C)^/^ + (1 ~ C)^^^]/2 being a spin-scaling 
factor. The correlation enhancement factor is 



F, 



1 + bisl + b24 
1 + bssl + bisi 



(21) 



with 61 = 1.01453936, 62 = 0.3255243, 63 = 0.941597104, 
and 64 = 0.587664306. Eq. ^ is a simple Fade approx- 
imation that recovers the RPA behavior of the uniform 
electron gas when Sc — 0- AH the parameters were found 
by the fitting procedure, and not by constraints on the 
integrated correlation energy (which would suggeslii^ an 
exponent of 7/6 and the appearance of (j) in the denom- 
inator of Eq. (j20p . and a quadratic term in the small- 
gradient expansion of Eq. jnH)). The irrelevance of some 
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FIG. 4: (ec - eJ^-^'*-^^^)/ec of the spin-unpolarized (C = 0) 
Airy gas model versus Sc (see Eq. (|20|l 'l for numerical RPA 
and our fit ARPA of Eq. (|19[l , for several slopes of the effective 
potential {F = 0.1, 0.5, and 1). Note that the numerical RPA 
has errors of order 2% in the region of small reduced gradient 

Sc. 
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FIG. 6: Enhancement factor F^c (see Eq. (|24|) for the spin- 
unpolarized case ((^ = 0), as a function of the reduced gradient 
s for several values of Va (ts = 0, 1, 2, 5, 10, and 20). The 
thin lines represent the ARPA+ enhancement factor whereas 
the thick lines are the PBEsol enhancement factor for rs = 
and Ts = 20 respectively. The LSDA is Fxc(ts, ^ = Oi * = 0). 



RPA fit (ARPA) for F=0.1, 0.5 and 1 




FIG. 5: (ec - e, 



LSDA-RPA 



)/ec of the fully-spin-polarized {C, = 



1) Airy gas model versus Sc for numerical RPA and our fit 
ARPA of Eq. (|19|l . for several slopes of the effective potential 
{F = 0.1, 0.5, and 1). Note that the numerical RPA has 
errors of order 2% in the region of small reduced gradient Sc. 



standard constraints may be related to the absence^ of 
a second-order gradient expansion for the conventional 
correlation energy density. Given F , e^^^ is a function 
of z, and Sc is a monotonic (hence invertible) function of 
z, so ef^^ can be expressed as a function of Sc- Since 
there is a one-to-one correspondence between the 
and our e^^^^, we can do the fitting. The fitting was 
done for Sc between and 20. 

In Figs. [3]and[n]we show (ec - e^^^^^^^^)/ec versus 
Sc for the spin-unpolarized Airy gas model (C = 0) and 
fully-spin-polarized Airy gas model (C = 1) respectively, 
for the slopes of the the effective potential used in Figs. [1] 
and[5](F = 0.1, 0.5, and 1). We note that our numerical 



calculation is accurate for Sc >^ 0.3, see Ref.'^^. We see 
in both figures that the numerical RPA correlation energy 
density does not depend much on the slope value F when 
they are plotted against Sc, motivating our definition of 
Sc in Eq. (|20|) and making the fit of the RPA correlation 
energy per particle independent of the F valued (see Eqs. 
dUl) and dni)). For Sc < 0.5 the ARPA of Eq. ^ is 
close to exact even if it does not match well the detailed 
exact behavior, as it does in the region 0.5 < Sc < 10. 
Overall we consider 



ARPA 



ARPA 



(22) 



an xc GGA functional that fits very well the Airy gas 
RPA xc energy density. Thus making the RPA -I- short- 
range correction (see Eq. to ARPA GGA, we pro- 
pose the following GGA xc functional (ARPA-I- GGA) 
constructed from the Airy gas 



^ARPA+ 



ARPA 



PBE 



PBE-RPA 



(23) 



The nonlocality of a GGA is displayed by the enhance- 
ment factor'''^''*'* 



GGA 



MGA 



uni f 



in) 



(24) 



gumf ^^-^ being the exchange energy per particle of a spin- 
unpolarized uniform electron gas. For a spin-unpolarized 
system in the high-density limit (rg — > 0), the exchange 
energy is dominant and Eq. defines the exchange 



enhancement factor = e^<^^(n, Vn)/e^™-''(n). 

Figs. [6] and [7] show the enhancement factor of ARPA-I- 
compared to PBEsol as a function of the reduced gradient 
s, for several values of r^, in the spin-unpolarized case 
and the fully-spin-polarized case, respectively. In both 
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FIG. 7: Enhancement factor Fxc (see Eq. (|2i|l for the fully- 
spin-polarized case = 1), as a function of the reduced gra- 
dient s for several values of Va (ts = 1, 2, 5, 10, and 20). The 
thin lines represent the ARPA-I- enhancement factor whereas 
the thick lines are the PBEsol enhancement factor for Ts = \ 
and Ts ~ 20 respectively. The LSD A is Fxdj's, = 1, s = 0). 
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FIG. 9: Comparison of F^J^^^+irsX = (shown with 

thin lines) and -/^^/^"^(r^, C = l,s) of Ref.ii (shown with 
thick lines) for several values of rs {ra = 1, 2, 5, and 20). 



IV. TESTS OF THE ARPA+ GGA XC ENERGY 
FUNCTIONAL 
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FIG. 8: Comparison of F^J^^^+{ra,C = 0, s) (shown with 
thin lines) and F^J''^'^^^ {rs, C, — 0, s) (shown with thick lines) 
for several values of rs (rs — 0, 2, 5, and 20). 



figures, the ARPA-I- and PBEsol enhancement factors 
agree well at small gradients (for s < 0.5), but for s >> 
0.5 ARPA-t- shows more exchange-correlation nonlocality 
than PBEsol. 

Figs. [5] and [5] show a comparison between the ARPA-I- 
GGA and AM05 GGA enhancement factors, for the spin- 
unpolarized and fully spin-polarized cases. Up to s = 0.5, 
F^/^^'^+irsXiS) and F^*^°^(rs, s) agree very well. 
For s > 0.5, F^^^'^+(rs, C, s) shows slightly more non- 
locality than i^^/^°^(r^,C,s), and, even if this difference 
is small, it has noticeable effects for the lattice constants 
of bulk solids. Overall, our ARPA-I- confirms the AM05 
construction for correlation. 



In this section we test our functionals for jellium sur- 
faces, atoms, molecules, and bulk solids. The calculations 
use the spin-scaling relation of Eq. . 



A. Jellium surfaces 

In Fig. [TO] we show e^f^ given by Eq. 1^, e^^^"^ 
given by Eq. dH]), and efi^^"^^^ of Rcf.^*, for two 
thick jellium slabs of bulk parameters rs — 2.07 and rs — 
4. We use accurate LSDA orbitals and densities as in 
Refs,S6d5i46. ARPA fits well the exact RPA until s w 20, 
showing that the Airy gas and the jellium surfaces are 
very close related, as expected. 

In Table [J we report the ARPA and ARPA-I- jellium 
surface exchange and xc energies. The a^'^^^^ are close 
to but worse than <7^^^ ■ The cr^^^^ are between cr^^^ 
and (7x^^^^^^ for rs <^ 3, but lower than both others 
for Tg >~ 4. The <Jxc^^^^ reasonably close to cr^^^^ 
and cr^^^ (see Ref4^), but are surprisingly lower and 
less accurate than a';:^^^. 



B. Spherical atoms 

In Table |TT] we calculate the ARPA-I- exchange and 
correlation energies of several atoms and ions. We use 
spin-restricted analytic Hartree-Fock orbitals'^^ and den- 
sities. (The difference between Hartree-Fock orbitals and 
Kohn-Sham orbitals is small for atoms.) For every atom 
and ion of Table |TT1 ARPA-H GGA improves the LSDA 
results, but it is still a poor approximation in comparison 
with GGA's constructed for atoms and molecules, such 
as PBE GGA^ii^. 
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FIG. 10: RPA exchange-correlation energy (hartree) per par- 
ticle e^c at position z versus z/Xf, at surlaces ol two jellium 
slabs. The bulk parameters are Vs = 2.07 and Vs = 4. Both 
jellium slabs have the width d = 3.2Af. The edges of the 
positive background are sd, z — 0. The differences at large 
z, emphasized here by plotting e^c instead of ne^c, are not 
important for the surface energy. 

In Table III we show the xc contribution to the valence- 
shell removal energy (a quantity that can be accurately 
measured experimentally^) of three atoms (Li, Be, and 
Ne). We observe that the ARPA+ systematically im- 
proves the LSDA results, competing in accuracy with the 
PBE GGA. 



C. Atomization energies of molecules 

The AE6 test set^^ of atomization energies of molecules 
has only six molecules (SiH4, SiO, S2, C3H4, C2H2O2, 
and C4H8) and was constructed to reproduce the errors 
of density functionals for larger molecular sets, providing 
a quick but representative evaluation of the accuracy of 
density functionals for molecules. In Table IIVI we show 
the errors (in kcal/mol) of the AE6 atomization ener- 
gies for ARPA-f GGA, ARPA GGA, PBE GGA, PBEsol 
GGA, and AMOS GGA. The errors given by ARPA+ 
GGA and ARPA GGA are practically the same, in ac- 
cord with the work of RefJ^, and show that the RPA+ 
short-range correction does not have an important effect 
on the atomization energies of molecules. Although our 
GGA short-range correction to RPA is important for to- 
tal energies, it tends to cancel out of energy differences 
for processes in which the electron number remains un- 
changed (as in Tables U and HV] but not Tables [TTl and 
The accuracy of the ARPA+ for the AEG test is close to 
that of PBEsol, with both reducing the LSDA error by 
by more than a factor of two. 

While our ARPA overbinds molecules (and this 
over binding is only slightly reduced in ARPA+), the full 
RPA apparently underbinds molecules^. Thus, even 
at the RPA level, the Airy gas xc energy density does 



TABLE I: Jellium surface exchange and exchange-correlation 
energies (erg/cm^) for LSDA, PBE, and ARPA-f in and be- 
yond the random phase approximation. We also show the jel- 
lium surface exchange and exchange-correlation energies be- 
yond RPA, for LAG GGA, AMOS GGA, PBEsol GGA, and 
TPSS meta-GGA of Ref.ii. The exact values of (7^'"=' and 
(J^^'^^ are from Ref.— , and the fixed-node difussion Monte 
Carlo (DMC) a^J^"^ values are interpolations and extrapola- 
tions of the estimates of Ref42i (see Table II of Ref. -'^ ) . To in- 
terpolate or extrapolate we recommend Eq. (15) of Ref.—. 
(Ihartree/bohr^ = 1.557 x lO^erg/cm^.) 
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not seem to transfer very accurately to molecules: much 
better atomization energies are predicted by standard 
functionals hke the PBE GGA^^ or the TPSS meta- 
GGAii. GGA overbinding of molecules typically goes 
together with GGA underestimation of the magnitude 
of the exchange-correlation energy of an atom, which we 
found for LSDA and ARPA+ but not so much for PBE 
in Table HH 



D. Equilibrium lattice constants of solids 

In Table|V]we test the ARPA+ GGA for a simple metal 
(Na), a semiconductor (Si), a transition metal (Cu), and 
an ionic solid (NaCl). The ARPA+ GGA lattice con- 
stants are longer than the PBEsol ones, but shorter than 
the PBE values, except for NaCl where ARPA -I- is close 
to PBE. These trends are plausible from the enhance- 
ment factors plotted in Figs. [6] and [71 and the maximum 
s values reported in Ref.'^'*. These calculations also sug- 



8 



TABLE 11: Exchange and correlation energies (in hartrees) of 
several spherical atoms and ions with spin-restricted Hartree- 
Fock orbitals and densitiesSi. Exact correlation energies are 
from Ref.— . PBE GGA, not shown in the table, has the mean 
absolute errors (m.a.e.): 0.0476 for exchange and 0.01563 for 
correlation. (See also Table V of Ref.=.) 





pLSDA 


T-,ARPA+ 




pLSDA 


j^ARPA+ 


^exact 


H 


-0.268 


-0.280 


-0.313 


-0.0222 


-0.0199 





He 


-0.884 


-0.925 


-1.026 


-0.1125 


-0.1030 


-0.0420 


Li+ 


-1.421 


-1.486 


-1.652 


-0.1346 


-0.1233 


-0.0435 


Be2+ 


-1.957 


-2.047 


-2.277 


-0.1504 


-0.1378 


-0.0443 


Li 


-1.538 


-1.603 


-1.781 


-0.1508 


-0.1378 


-0.0453 


Be+ 


-2.168 


-2.261 


-2.507 


-0.1727 


-0.1578 


-0.0474 


Be 


-2.312 


-2.408 


-2.667 


-0.2240 


-0.2058 


-0.0943 


B+ 


-3.036 


-3.157 


-3.492 


-0.2520 


-0.2317 


-0.1113 


Nc6+ 


-6.634 


-6.886 


-7.594 


-0.3336 


-0.3069 


-0.1799 


N 


-5.893 


-6.047 


-6.596 


-0.4273 


-0.4016 


-0.1883 


Ne 


-11.033 


-11.220 


-12.109 


-0.7428 


-0.7084 


-0.3905 


Ar 


-27.863 


-28.118 


-30.190 


-1.4242 


-1.3723 


-0.7222 


m.a.e. 


0.600 


0.481 




0.1865 


0.1664 





TABLE IIL Change in xc energy (hartree) of an atom due to 
removal of a sheU of valence electrons(A£;i:c = E^l""" - E^°"). 
The calculation is based on the exchange and correlation en- 
ergies listed in Table |ll] of this work and in Table VI of Ref.—. 











/\ jpexact 


Li Li+ 
Be ^ Be+2 
Ne Ne+"5 


-0.133 
-0.429 
-4.808 


-0.132 
-0.430 
-4.737 


-0.138 
-0.438 
-4.793 


-0.131 
-0.440 
-4.726 



TABLE IV: The errors (kcal/mole) of the atomization en- 
ergies of the AE6 set of molecules. We use the 6 — 311 + 
G(3d/, 2p) basis set in the GaussianOS code. The AMOS at- 
omization energies of the AE6 set of molecules were calculated 
in Ref.—, using the spin-polarized version of AM05 given in 
Reffil. The LSDA mean error (ME) is 77.3 kcal/mole and 
its mean absolute error (MAE) is 77.3 kcal/mole^^ The 
TPSS meta-GGA of Ref.-^ gives ME=4.2 kcal/mole, and 
MAE=6.0 kcal/mole. The AE6 mean atomization energy is 
517 kcal/mole. (1 hartree = 627.5 kcal/mole.) (For ARPA+ 
and ARPA, we used PBEsol densities.) 





PBE 


ARPA-(- 


ARPA 


PBEsol 


AM05 


SiH4 


-9.2 


10.1 


9.9 


1.3 


7.6 


SiO 


3.6 


11.2 


12.3 


12.9 


13.5 


S2 


13.1 


18.4 


19.2 


21.9 


21.6 


C3H4 


16.4 


46.0 


50.6 


45.1 


48.1 


C2H2O2 


31.8 


60.1 


65.7 


64.7 


66.6 


C4H8 


18.7 


70.6 


78.7 


69.6 


75.0 


ME 


12.4 


36.1 


39.4 


35.9 


38.7 


MAE 


15.5 


36.1 


39.4 


35.9 


38.7 



gest that the correct second-order gradient expansion for 
exchange^, employed in the construction of the PBEsol 
GGA, is the most promising path toward an accurate and 
noncmpirical GGA for soUds. 

TABLE V: Lattice constants (in A) calculated with the Gaus- 
sian03 code as in Ref.— and compared to experimental values 
corrected to the static-lattice limit^!^-. (For ARPA+, we 
used PBEsol densities.) 



Solid 


LSDA 


PBE 


PBEsol 


ARPA-I- 


Exp or. 


Na 


4.049 


4.199 


4.159 


4.207 


4.210 


Si 


5.410 


5.479 


5.442 


5.470 


5.423 


Cu 


3.530 


3.635 


3.578 


3.605 


3.596 


NaCl 


5.471 


5.696 


5.611 


5.716 


5.580 


ME 


-0.087 


0.050 


-0.005 


0.045 




MAE 


0.087 


0.056 


0.030 


0.049 





The GaussianOS code that we use gives lattice con- 
stants that are on average a little too long^''. The LSDA 
lattice constants calculated with the more-accurate 
WIEN2K code ar(^: Na 4.047, Si 5.407, Cu 3.522, and 
NaCl 5.465. Thus, extensive and more accurate lat- 
tice constants calculations need to be performed for our 
ARPA+. 



V. CONCLUSIONS 

In this paper we construct the RPA correlation energy 
density of the Airy gas, using an accurate Airy gas model 
that has only 19 occupied orbitals. This approch can be 
generalized to other physical systems, such as a more so- 
phisticated edge electron gas that can include curvature 
corrections (arising from nonlinearity of Veff{z)). 

We have constructed the ARPA GGA that accurately 
fits the RPA xc energy density of the Airy gas, and 
we have corrected its short-range part in the framework 
of the RPA-I- 15 approach, developing the ARPA -I- GGA 
entirely without empiricism. Because of the delocaliza- 
tion of the electrons in the Airy gas, our ARPA-f GGA 
has nearly the correct Airy-gas correlation energy. Via 
our Figs. 8 and 9, our ARPA-f confirms the AM05 
hypothesisi^ for the correlation functional compatible 
with Airy-gas GGA exchange^'^'^'*. 

By testing the ARPA-I- GGA for jellium surfaces, 
atoms, molecules, and bulk solids, we observe that the 
xc energy density of the Airy gas can be transferred suc- 
cessfully to a very similar system such as the jellium sur- 
face, but less successfully to a very different system like a 
bulk solid, an atom, or a molecule. However, the ARPA-I- 
GGA mildly improves the LSDA results for lattice con- 
stants and atomization energies, without much worsening 
the already-good surface exchange-correlation energies. 

We would have liked to replace the RPA-I- method 
by the more sophisticated inhomogeneous Singwi-Tosi- 
Land-Sjolander (ISTLS)i^i^, but were not able to 
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achieve sufficiently accurate numerical results for the cor- 
relation energy densities thereof. The future use of ISTLS 
could refine our input, and provide an energy density 
(not just an integrated energy) for the short-range cor- 
rection to RPA. Other possible future refinements could 
include the use of different reference systems for the 
bulk and surface of a solidliiAS,, replacing the Airy gas 
by a more sophisticated example of the edge electron 
gas, or replacing the GGA functional form by the meta- 
GGAil. We suspect^ii^ that the meta-GGA form is 
needed to achieve simultaneous high accuracy for atoms, 
molecules, and solids near equilibrium. In fact the TPSS 
meta-GG Ai^"'?^ is already close to being such a general- 
purpose semilocal functional, and a revised TPSS^ with 
improved lattice constants may be even closer. 

We note however that there are two formally unsatis- 
factory aspects of using the exchange-correlation energy 
density of a nonuniform system as a reference for the 
construction of density functionals: (1) Except in the 
uniform electron gas, the energy density is neither ob- 
servable nor unique, since any function integrating to 
zero can be added to it with no physical consequence. 
Here, as in RefsJ^fi^^r^, and^, we have chosen the 
conventional^ gauge for the energy density, but other 
choices should be explored. (2) While the integrated ex- 
change energy for a slowly- varying density is expressible 



in terms of the GGA ingredients n and Vn, the con- 
ventional exchange energy density in this limit is not so 
expressible, having a Laplacian term V^n^/^ which inte- 
grates to zero but has a divergent coefficien t ""^^i^^ . As a re- 
sult, the Airy-gas GGA cannot predict accurate exchange 
energies for slowly- varying electron densities (e.g., the jel- 
lium surface exchange energy), while more standardly- 
constructed GGA's like PBEsol can do sc^^ (our Table 
I). The Airy-gas GGA can at best work for the jellium 
surface by error cancellation between exchange and corre- 
lation, which is possible for typical valence-electron den- 
sities but not under uniform density scaling to the high- 
density limit where exchange dominates. 

The GGA constructed here has no clear practical ad- 
vantage over already-published ones. Our purpose is not 
to advocate its use, but to show what is obtained from 
the Airy-gas reference system within a consistent imple- 
mentation for correlation as well as exchange. 
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